## World War 2 Manuscript Dataset ##
## Prepared by Kevin Fahey ##
## And Doug Atkinson ##
## Compiled with R 4.2.1 ##

##### Preamble #####
rm(list=ls())
set.seed(20120630)

# Load packages
library(foreign)
library(sandwich)
library(stargazer); library(xtable)
library(tidyverse); library(broom); library(dotwhisker)
library(fixest)
library(haven)
library(plm)
library(pglm)
library(stringr)
library(maptools); library(usmap)
library(foreach)
library(zoo)
library(multiwayvcov); library(lmtest); library(lfe)
library(psych); library(knitr)
library(gridExtra)
options(scipen=7)
options(digits=4)

# Set working directory here!


## Load in datasets ##
dat <- read.csv("dat.csv")
dat.analysis <- read.csv("datanalysis.csv")
final.statedat <- read.csv("finalstatedat.csv")
statedat <- read.csv("statedat.csv")


## Create annual datasets ##
dat42 <- dat.analysis %>% filter(year == 42)
dat43 <- dat.analysis %>% filter(year == 43)
dat44 <- dat.analysis %>% filter(year == 44)
dat45 <- dat.analysis %>% filter(year == 45)


##### Table 1 -- Descriptive Statistics #####
dat_table1 <- dat.analysis %>%
  dplyr::select(propelig, propdraft, male1554.elig, draft_melig,
                swingcounty, splitswing, pickupswing, swingcongavg, 
                GOPAvg, urb940, totpop40, cropval, 
                m1554.interpolated, pctblack, farmprop) %>%
  mutate(urb940 = urb940/totpop40,
         cropval = cropval/totpop40,
         totpop40 = totpop40/1000) %>%
  rename(Enlistees_Population = propelig,
         Draftees_Population = propdraft,
         Enlistees_Male_Eligible = male1554.elig,
         Draftees_Male_Eligible = draft_melig,
         Swing_County_50_55 = swingcounty,
         Swing_County_47.5_53.5 = splitswing,
         Swing_County_45_50 = pickupswing,
         Swing_County_Congress = swingcongavg,
         Pct_GOP = GOPAvg,
         Pct_Urban_1940 = urb940,
         Total_Population = totpop40,
         Per_Capita_Crop_Value = cropval,
         Pct_Male_Eligible = m1554.interpolated,
         Pct_Blank = pctblack,
         Farms_Population = farmprop) %>%
  psych::describe(., skew = F, ranges = T) %>%
  dplyr::select(-vars, -range, -se) %>%
  dplyr::rename(Obs = n,
                Average = mean,
                Standard_Deviation = sd,
                Min = min,
                Max = max) %>%
  dplyr::mutate(Obs = round(Obs, 0),
                Average = round(Average, 3),
                Standard_Deviation = round(Standard_Deviation, 3),
                Min = round(Min, 3),
                Max = round(Max, 3)) %>%
  dplyr::select(Average, Standard_Deviation, Min, Max, Obs) %>%
  knitr::kable(align = "c", format = "latex",
               caption = "Summary Statistics, covariates. Each observation represents a single county for a single year."); dat_table1

##### Table 2 -- Main Models #####

# original/baseline models
main_all <- feols(propelig ~ swingcounty
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop
            + yr42 + yr43 + yr44 + yr45,
            data = dat.analysis,
            cluster = "fips")
main_elec <- feols(propelig ~ swingcounty + electionyear
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat.analysis,
            cluster = "fips")
main_inter <- feols(propelig ~ swingcounty * electionyear
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat.analysis,
            cluster = "fips")
main_42 <- feols(propelig ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop, 
            data = dat42,
            se = "hetero")
main_43 <- feols(propelig ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop, 
            data = dat43,
            se = "hetero")
main_44 <- feols(propelig ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop, 
            data = dat44,
            se = "hetero")
main_45 <- feols(propelig ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop, 
            data = dat45,
            se = "hetero")
main_am1 <- feols(propelig ~ splitswing
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + yr42 + yr43 + yr44 + yr45, 
             data = dat.analysis,
             cluster = "fips")
main_am2 <- feols(propelig ~ pickupswing
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + yr42 + yr43 + yr44 + yr45, 
             data = dat.analysis,
             cluster = "fips")
main_am3 <- feols(propelig ~ swing40 
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + yr42 + yr43 + yr44 + yr45, 
             data = dat.analysis,
             cluster = "fips")

# draftees only
draft_all <- feols(propdraft ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45,
                  data = dat.analysis,
                  cluster = "fips")
draft_elec <- feols(propdraft ~ swingcounty + electionyear
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop,
                   data = dat.analysis,
                   cluster = "fips")
draft_inter <- feols(propdraft ~ swingcounty * electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
draft_42 <- feols(propdraft ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat42,
                 se = "hetero")
draft_43 <- feols(propdraft ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat43,
                 se = "hetero")
draft_44 <- feols(propdraft ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat44,
                 se = "hetero")
draft_45 <- feols(propdraft ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat45,
                 se = "hetero")
draft_am1 <- feols(propdraft ~ splitswing
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45, 
                  data = dat.analysis,
                  cluster = "fips")
draft_am2 <- feols(propdraft ~ pickupswing
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45, 
                  data = dat.analysis,
                  cluster = "fips")
draft_am3 <- feols(propdraft ~ swing40 
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45, 
                  data = dat.analysis,
                  cluster = "fips")

# male-eligible interpolated DV
melig_all <- feols(male1554.elig ~ swingcounty
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45,
                   data = dat.analysis,
                   cluster = "fips")
melig_elec <- feols(male1554.elig ~ swingcounty + electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
melig_inter <- feols(male1554.elig ~ swingcounty * electionyear
                     + GOPAvg + pcturban + totpop.thou
                     + avgcropvalue + pctmale.elig2
                     + pctblack + farmprop,
                     data = dat.analysis,
                     cluster = "fips")
melig_42 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat42,
                  se = "hetero")
melig_43 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat43,
                  se = "hetero")
melig_44 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat44,
                  se = "hetero")
melig_45 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat45,
                  se = "hetero")
melig_am1 <- feols(male1554.elig ~ splitswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
melig_am2 <- feols(male1554.elig ~ pickupswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
melig_am3 <- feols(male1554.elig ~ swing40 
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")

# male-eligible interpolated DV, draftees only
delig_all <- feols(draft_melig ~ swingcounty
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45,
                   data = dat.analysis,
                   cluster = "fips")
delig_elec <- feols(draft_melig ~ swingcounty + electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
delig_inter <- feols(draft_melig ~ swingcounty * electionyear
                     + GOPAvg + pcturban + totpop.thou
                     + avgcropvalue + pctmale.elig2
                     + pctblack + farmprop,
                     data = dat.analysis,
                     cluster = "fips")
delig_42 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat42,
                  se = "hetero")
delig_43 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat43,
                  se = "hetero")
delig_44 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat44,
                  se = "hetero")
delig_45 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat45,
                  se = "hetero")
delig_am1 <- feols(draft_melig ~ splitswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
delig_am2 <- feols(draft_melig ~ pickupswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
delig_am3 <- feols(draft_melig ~ swing40 
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")


## Produce coefficient plot of all models
coefplot_1 <- rbind(tidy(main_all) %>%
                       filter(term == "swingcounty") %>%
                       mutate(model = "Baseline"),
                     tidy(main_elec) %>%
                       filter(term == "swingcounty") %>%
                       mutate(model = "Election Dummy"),
                     tidy(main_inter) %>%
                       filter(term == "swingcounty:electionyear") %>%
                       mutate(model = "Election Interaction"),
                     tidy(main_42) %>%
                       filter(term == "swingcounty") %>%
                       mutate(model = "1942 Only"),
                     tidy(main_43) %>%
                       filter(term == "swingcounty") %>%
                       mutate(model = "1943 Only"),
                     tidy(main_44) %>%
                       filter(term == "swingcounty") %>%
                       mutate(model = "1944 Only"),
                    tidy(main_45) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1945 Only"),
                     tidy(main_am1) %>%
                       filter(term == "splitswing") %>%
                       mutate(model = "Swing 47.5-52.5"),
                     tidy(main_am2) %>%
                       filter(term == "pickupswing") %>%
                       mutate(model = "Swing 45-50")) %>%
  mutate(term = model)

coefplot_2 <- rbind(tidy(draft_all) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "Baseline"),
                    tidy(draft_elec) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "Election Dummy"),
                    tidy(draft_inter) %>%
                      filter(term == "swingcounty:electionyear") %>%
                      mutate(model = "Election Interaction"),
                    tidy(draft_42) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1942 Only"),
                    tidy(draft_43) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1943 Only"),
                    tidy(draft_44) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1944 Only"),
                    tidy(draft_45) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1945 Only"),
                    tidy(draft_am1) %>%
                      filter(term == "splitswing") %>%
                      mutate(model = "Swing 47.5-52.5"),
                    tidy(draft_am2) %>%
                      filter(term == "pickupswing") %>%
                      mutate(model = "Swing 45-50")
) %>%
  mutate(term = model)

coefplot_3 <- rbind(tidy(melig_all) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "Baseline"),
                    tidy(melig_elec) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "Election Dummy"),
                    tidy(melig_inter) %>%
                      filter(term == "swingcounty:electionyear") %>%
                      mutate(model = "Election Interaction"),
                    tidy(melig_42) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1942 Only"),
                    tidy(melig_43) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1943 Only"),
                    tidy(melig_44) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1944 Only"),
                    tidy(melig_45) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1945 Only"),
                    tidy(melig_am1) %>%
                      filter(term == "splitswing") %>%
                      mutate(model = "Swing 47.5-52.5"),
                    tidy(melig_am2) %>%
                      filter(term == "pickupswing") %>%
                      mutate(model = "Swing 45-50")
) %>%
  mutate(term = model)

coefplot_4 <- rbind(tidy(delig_all) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "Baseline"),
                    tidy(delig_elec) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "Election Dummy"),
                    tidy(delig_inter) %>%
                      filter(term == "swingcounty:electionyear") %>%
                      mutate(model = "Election Interaction"),
                    tidy(delig_42) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1942 Only"),
                    tidy(delig_43) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1943 Only"),
                    tidy(delig_44) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1944 Only"),
                    tidy(delig_45) %>%
                      filter(term == "swingcounty") %>%
                      mutate(model = "1945 Data Only"),
                    tidy(delig_am1) %>%
                      filter(term == "splitswing") %>%
                      mutate(model = "Swing 47.5-52.5"),
                    tidy(delig_am2) %>%
                      filter(term == "pickupswing") %>%
                      mutate(model = "Swing 45-50")
) %>%
  mutate(term = model)


# output as plots
library(dotwhisker)
models_output_1 <- dwplot(coefplot_1, dodge_size = 0.5, ci = 0.95,
                        vline = geom_vline(xintercept = 0.0,
                                           colour = "black",
                                           linetype = 2)) +
  theme_bw() + xlab("Enlistees/Population") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "none") + 
  scale_colour_grey(start = .3, end = .4) +
  guides(
    colour = guide_legend("Model"),
    linetype = guide_legend("Model")
  )
png("models_output_1.png", res = 100)
print(models_output_1)
dev.off()

models_output_2 <- dwplot(coefplot_2, dodge_size = 0.5, ci = 0.95,
                          vline = geom_vline(xintercept = 0.0,
                                             colour = "black",
                                             linetype = 2)) +
  theme_bw() + xlab("Draftees/Population") + ylab("") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "none") + 
  scale_y_discrete(position = "right") + 
  scale_colour_grey(start = .3, end = .4) +
  guides(
    colour = guide_legend("Model"),
    linetype = guide_legend("Model")
  )
png("models_output_2.png", res = 100)
print(models_output_2)
dev.off()

models_output_3 <- dwplot(coefplot_3, dodge_size = 0.5, ci = 0.95,
                          vline = geom_vline(xintercept = 0.0,
                                             colour = "black",
                                             linetype = 2)) +
  theme_bw() + xlab("Enlistees/Male Eligible Population") + ylab("") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "none") + 
  scale_colour_grey(start = .3, end = .4) +
  guides(
    colour = guide_legend("Model"),
    linetype = guide_legend("Model")
  )
png("models_output_3.png", res = 100)
print(models_output_3)
dev.off()


models_output_4 <- dwplot(coefplot_4, dodge_size = 0.5, ci = 0.95,
                          vline = geom_vline(xintercept = 0.0,
                                             colour = "black",
                                             linetype = 2)) +
  theme_bw() + xlab("Draftees/Male Eligible Population") + ylab("") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "none") + 
  scale_y_discrete(position = "right") + 
  scale_colour_grey(start = .3, end = .4) +
  guides(
    colour = guide_legend("Model"),
    linetype = guide_legend("Model")
  )
png("models_output_4.png", res = 100)
print(models_output_4)
dev.off()


# E-Tables
etable(m1, m2, m3, m42, m43, m44, m45, 
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Impact of county partisanship on U.S. Army enlistment rates in World War 2. Unit-clustered standard errors listed in parentheses.",
       label = "reg1",
       fitstat = c("n", "ar2", "f"))
etable(am1, am2, am3,
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Impact of alternate specifications of swing-county status on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered fixed effects employed in all models.",
       label = "robust1",
       fitstat = c("n", "ar2", "f"))
etable(statem1, statem2, statem3, statem4,
       cluster = "statefips", se = "hetero", tex = TRUE,
       drop = "statefips", digits = "r3", digits.stats = "r3",
       title = "Impact of state presidential vote share on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered standard employed in all models.",
       label = "stateagg",
       fitstat = c("n", "ar2", "f"))


##### Aggregate to State #####
statem1 <- feols(enlistpop ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem2 <- feols(enlistpop ~ swing36
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem3 <- feols(enlistpop ~ swing40
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem4 <- feols(enlistpop ~ swingavg
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")

statem5 <- feols(enlistmale ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem6 <- feols(enlistmale ~ swing36
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem7 <- feols(enlistmale ~ swing40
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem8 <- feols(enlistmale ~ swingavg
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")

statem9 <- feols(draftpop ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem10 <- feols(draftpop ~ swing36
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem11 <- feols(draftpop ~ swing40
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem12 <- feols(draftpop ~ swingavg
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")

statem13 <- feols(draftmale ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem14 <- feols(draftmale ~ swing36
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem15 <- feols(draftmale ~ swing40
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem16 <- feols(draftmale ~ swingavg
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")


# plot
tidystate <- rbind(tidy(statem1) %>%
                     filter(term == "swing32") %>% 
                     mutate(model = "1932 Election"),
  tidy(statem2) %>%
    filter(term == "swing36") %>%
    mutate(model = "1936 Election"),
  tidy(statem3) %>%
    filter(term == "swing40") %>%
    mutate(model = "1940 Election"),
  tidy(statem4) %>%
    filter(term == "swingavg") %>%
    mutate(model = "Average of Elections"),
  tidy(statem5) %>%
    filter(term == "swing32") %>% 
    mutate(model = "1932 Election"),
  tidy(statem6) %>%
    filter(term == "swing36") %>%
    mutate(model = "1936 Election"),
  tidy(statem7) %>%
    filter(term == "swing40") %>%
    mutate(model = "1940 Election"),
  tidy(statem8) %>%
    filter(term == "swingavg") %>%
    mutate(model = "Average of Elections"),
  tidy(statem9) %>%
    filter(term == "swing32") %>% 
    mutate(model = "1932 Election"),
  tidy(statem10) %>%
    filter(term == "swing36") %>%
    mutate(model = "1936 Election"),
  tidy(statem11) %>%
    filter(term == "swing40") %>%
    mutate(model = "1940 Election"),
  tidy(statem12) %>%
    filter(term == "swingavg") %>%
    mutate(model = "Average of Elections"),
  tidy(statem13) %>%
    filter(term == "swing32") %>% 
    mutate(model = "1932 Election"),
  tidy(statem14) %>%
    filter(term == "swing36") %>%
    mutate(model = "1936 Election"),
  tidy(statem15) %>%
    filter(term == "swing40") %>%
    mutate(model = "1940 Election"),
  tidy(statem16) %>%
    filter(term == "swingavg") %>%
    mutate(model = "Average of Elections")
  ) %>%
  mutate(term = model,
         model = c(rep("Enlist/Pop", 4),
                   rep("Enlist/M. Elig", 4),
                   rep("Draft/Pop", 4),
                   rep("Draft/M. Elig", 4)))


# coefplot output
stateaggregate <- dwplot(tidystate, dodge_size = 0.5, ci = 0.95,
                         dot_args = list(aes(shape = model)),
                         whisker_args = list(aes(linetype = model)),
                         vline = geom_vline(xintercept = 0,
                                            colour = "grey60",
                                            linetype = 2)) +
  theme_bw() + xlab("Coefficient Estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom") + 
 guides(colour = guide_legend(override.aes = list(shape = c(0, 1, 2, 3), 
                                                  linetype = c(0, 1, 2, 3)))) + 
 scale_shape_manual(name = "Model",
                      values = c(0, 1, 2, 3),
                    guide = "none") +
  scale_linetype_manual(name = "Model",
                        values = c(0, 1, 2, 3),
                        guide = "none") + 
  scale_colour_grey(start = .2, end = .25,
                    guide = "legend")
png("stateaggregate.png", res = 80)
print(stateaggregate)
dev.off()


##### Figure 4 -- Placebo Tests #####
## Five-pp bins from 30 to 70 percent Democrat ## 
dat.analysis$placebo1 <- ifelse(dat.analysis$DemAvg < 35 & dat.analysis$DemAvg >= 30, 1, 0)
dat.analysis$placebo2 <- ifelse(dat.analysis$DemAvg < 40 & dat.analysis$DemAvg >= 35, 1, 0)
dat.analysis$placebo3 <- ifelse(dat.analysis$DemAvg < 45 & dat.analysis$DemAvg >= 40, 1, 0)
dat.analysis$placebo4 <- ifelse(dat.analysis$DemAvg < 50 & dat.analysis$DemAvg >= 45, 1, 0)
dat.analysis$placebo5 <- ifelse(dat.analysis$DemAvg < 55 & dat.analysis$DemAvg >= 50, 1, 0)
dat.analysis$placebo6 <- ifelse(dat.analysis$DemAvg < 60 & dat.analysis$DemAvg >= 55, 1, 0)
dat.analysis$placebo7 <- ifelse(dat.analysis$DemAvg < 65 & dat.analysis$DemAvg >= 60, 1, 0)
dat.analysis$placebo8 <- ifelse(dat.analysis$DemAvg < 70 & dat.analysis$DemAvg >= 65, 1, 0)
dat.analysis$placebo9 <- ifelse(dat.analysis$DemAvg < 75 & dat.analysis$DemAvg >= 70, 1, 0)

## Models for output ##
placebo.m1 <- feols(propelig ~ placebo1
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m2 <- feols(propelig ~ placebo2
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m3 <- feols(propelig ~ placebo3
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m4 <- feols(propelig ~ placebo4
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m5 <- feols(propelig ~ placebo5
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m6 <- feols(propelig ~ placebo6
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m7 <- feols(propelig ~ placebo7
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m8 <- feols(propelig ~ placebo8
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m9 <- feols(propelig ~ placebo9
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")

## Assemble coefficients, ses, cis, and xlabs ##
out_placebo_a <- rbind(tidy(placebo.m1) %>%
               filter(term == "placebo1") %>%
               mutate(model = "30-35"),
             tidy(placebo.m2) %>%
               filter(term == "placebo2") %>%
               mutate(model = "35-40"),
             tidy(placebo.m3) %>%
               filter(term == "placebo3") %>%
               mutate(model = "40-45"),
             tidy(placebo.m4) %>%
               filter(term == "placebo4") %>%
               mutate(model = "45-50"),
             tidy(placebo.m5) %>%
               filter(term == "placebo5") %>%
               mutate(model = "50-55"),
             tidy(placebo.m6) %>%
               filter(term == "placebo6") %>%
               mutate(model = "50-60"),
             tidy(placebo.m7) %>%
               filter(term == "placebo7") %>%
               mutate(model = "60-65"),
             tidy(placebo.m8) %>%
               filter(term == "placebo8") %>%
               mutate(model = "65-70"),
             tidy(placebo.m9) %>%
               filter(term == "placebo9") %>%
               mutate(model = "70-75"))
coefs <- out_placebo_a$estimate
ses <- out_placebo_a$std.error
ci.lo <- coefs - (ses * 1.96)
ci.hi <- coefs + (ses * 1.96)
xlabs <-  out_placebo_a$model

## And plot ##
placebo.test <- plot(coefs, pch = 16, col = "black", lwd = 6,
                     main = "", ylab = "Coefficient, Placebo Test", ylim = c(-0.25, 0.4),
                     xlab = "Democratic Vote Share (Five Percentage Point Bins)",
                     xaxt="n")
abline(h = 0, col = "black", lwd = 2, lty = 3)
segments(x0 = 1, y0 = ci.lo[1], x1 = 1, y1 = ci.hi[1],
         col = "darkgrey", lwd = 3)
segments(x0 = 2, y0 = ci.lo[2], x1 = 2, y1 = ci.hi[2],
         col = "darkgrey", lwd = 3)
segments(x0 = 3, y0 = ci.lo[3], x1 = 3, y1 = ci.hi[3],
         col = "darkgrey", lwd = 3)
segments(x0 = 4, y0 = ci.lo[4], x1 = 4, y1 = ci.hi[4],
         col = "darkgrey", lwd = 3)
segments(x0 = 5, y0 = ci.lo[5], x1 = 5, y1 = ci.hi[5],
         col = "darkgrey", lwd = 3)
segments(x0 = 6, y0 = ci.lo[6], x1 = 6, y1 = ci.hi[6],
         col = "darkgrey", lwd = 3)
segments(x0 = 7, y0 = ci.lo[7], x1 = 7, y1 = ci.hi[7],
         col = "darkgrey", lwd = 3)
segments(x0 = 8, y0 = ci.lo[8], x1 = 8, y1 = ci.hi[8],
         col = "darkgrey", lwd = 3)
segments(x0 = 9, y0 = ci.lo[9], x1 = 9, y1 = ci.hi[9],
         col = "darkgrey", lwd = 3)
axis(1, seq(1, 9, 1), labels = xlabs)

# print #
dev.copy(png, 'placebotest.png')
dev.off()


## Now 2-pp bins from 40 to 60 percent ##
dat.analysis$placebo1 <- ifelse(dat.analysis$DemAvg < 42 & dat.analysis$DemAvg >= 40, 1, 0)
dat.analysis$placebo2 <- ifelse(dat.analysis$DemAvg < 44 & dat.analysis$DemAvg >= 42, 1, 0)
dat.analysis$placebo3 <- ifelse(dat.analysis$DemAvg < 46 & dat.analysis$DemAvg >= 44, 1, 0)
dat.analysis$placebo4 <- ifelse(dat.analysis$DemAvg < 48 & dat.analysis$DemAvg >= 46, 1, 0)
dat.analysis$placebo5 <- ifelse(dat.analysis$DemAvg < 50 & dat.analysis$DemAvg >= 48, 1, 0)
dat.analysis$placebo6 <- ifelse(dat.analysis$DemAvg < 52 & dat.analysis$DemAvg >= 50, 1, 0)
dat.analysis$placebo7 <- ifelse(dat.analysis$DemAvg < 54 & dat.analysis$DemAvg >= 52, 1, 0)
dat.analysis$placebo8 <- ifelse(dat.analysis$DemAvg < 56 & dat.analysis$DemAvg >= 54, 1, 0)
dat.analysis$placebo9 <- ifelse(dat.analysis$DemAvg < 58 & dat.analysis$DemAvg >= 56, 1, 0)
dat.analysis$placebo10 <- ifelse(dat.analysis$DemAvg < 60 & dat.analysis$DemAvg >= 58, 1, 0)

## Models for output ##
placebo.m1 <- feols(propelig ~ placebo1
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m2 <- feols(propelig ~ placebo2
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m3 <- feols(propelig ~ placebo3
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m4 <- feols(propelig ~ placebo4
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m5 <- feols(propelig ~ placebo5
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m6 <- feols(propelig ~ placebo6
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m7 <- feols(propelig ~ placebo7
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m8 <- feols(propelig ~ placebo8
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m9 <- feols(propelig ~ placebo9
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
placebo.m10 <- feols(propelig ~ placebo10
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")

## Assemble coefficients, ses, cis, and xlabs ##
out_placebo_b <- rbind(tidy(placebo.m1) %>%
                         filter(term == "placebo1") %>%
                         mutate(model = "40-42"),
                       tidy(placebo.m2) %>%
                         filter(term == "placebo2") %>%
                         mutate(model = "42-44"),
                       tidy(placebo.m3) %>%
                         filter(term == "placebo3") %>%
                         mutate(model = "44-46"),
                       tidy(placebo.m4) %>%
                         filter(term == "placebo4") %>%
                         mutate(model = "46-48"),
                       tidy(placebo.m5) %>%
                         filter(term == "placebo5") %>%
                         mutate(model = "48-50"),
                       tidy(placebo.m6) %>%
                         filter(term == "placebo6") %>%
                         mutate(model = "50-52"),
                       tidy(placebo.m7) %>%
                         filter(term == "placebo7") %>%
                         mutate(model = "52-54"),
                       tidy(placebo.m8) %>%
                         filter(term == "placebo8") %>%
                         mutate(model = "54-56"),
                       tidy(placebo.m9) %>%
                         filter(term == "placebo9") %>%
                         mutate(model = "56-58"),
                       tidy(placebo.m10) %>%
                         filter(term == "placebo10") %>%
                         mutate(model = "58-60"))
coefs.trunc <- out_placebo_b$estimate
ses.trunc <- out_placebo_b$std.error
ci.lo.trunc <- coefs.trunc - (ses.trunc * 1.96)
ci.hi.trunc <- coefs.trunc + (ses.trunc * 1.96)
xlabs.trunc <-  out_placebo_b$model

## And plot ##
placebo.test.trunc <- plot(coefs.trunc, pch = 16, col = "black", lwd = 6,
                           main = "", ylab = "Coefficient, Placebo Test", ylim = c(-0.2, 0.2),
                           xlab = "Democratic Vote Share (Two Percentage Point Bins)",
                           xaxt="n")
abline(h = 0, col = "black", lwd = 2, lty = 3)
segments(x0 = 1, y0 = ci.lo.trunc[1], x1 = 1, y1 = ci.hi.trunc[1],
         col = "darkgrey", lwd = 3)
segments(x0 = 2, y0 = ci.lo.trunc[2], x1 = 2, y1 = ci.hi.trunc[2],
         col = "darkgrey", lwd = 3)
segments(x0 = 3, y0 = ci.lo.trunc[3], x1 = 3, y1 = ci.hi.trunc[3],
         col = "darkgrey", lwd = 3)
segments(x0 = 4, y0 = ci.lo.trunc[4], x1 = 4, y1 = ci.hi.trunc[4],
         col = "darkgrey", lwd = 3)
segments(x0 = 5, y0 = ci.lo.trunc[5], x1 = 5, y1 = ci.hi.trunc[5],
         col = "darkgrey", lwd = 3)
segments(x0 = 6, y0 = ci.lo.trunc[6], x1 = 6, y1 = ci.hi.trunc[6],
         col = "darkgrey", lwd = 3)
segments(x0 = 7, y0 = ci.lo.trunc[7], x1 = 7, y1 = ci.hi.trunc[7],
         col = "darkgrey", lwd = 3)
segments(x0 = 8, y0 = ci.lo.trunc[8], x1 = 8, y1 = ci.hi.trunc[8],
         col = "darkgrey", lwd = 3)
segments(x0 = 9, y0 = ci.lo.trunc[9], x1 = 9, y1 = ci.hi.trunc[9],
         col = "darkgrey", lwd = 3)
segments(x0 = 10, y0 = ci.lo.trunc[10], x1 = 10, y1 = ci.hi.trunc[10],
         col = "darkgrey", lwd = 3)
axis(1, seq(1, 10, 1), labels = xlabs.trunc)

# print #
dev.copy(png, 'placebotesttrunc.png')
dev.off()


##### Figure 5 -- Congressional Results #####
mcavg <- lm(propelig ~ swingcongavg
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat)
se.mcavg <- sqrt(diag(vcovHC(mcavg)))

mc36a <- lm(propelig ~ swingcong36
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc36a <- sqrt(diag(vcovHC(mc36a)))

mc36b <- lm(propelig ~ splitcongswing36
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc36b <- sqrt(diag(vcovHC(mc36b)))

mc36c <- lm(propelig ~ pickupcongswing36
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc36c <- sqrt(diag(vcovHC(mc36c)))

mc38a <- lm(propelig ~ swingcong38
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc38a <- sqrt(diag(vcovHC(mc38a)))

mc38b <- lm(propelig ~ splitcongswing38
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc38b <- sqrt(diag(vcovHC(mc38b)))

mc38c <- lm(propelig ~ pickupcongswing38
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc38c <- sqrt(diag(vcovHC(mc38c)))

mc40a <- lm(propelig ~ swingcong40
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc40a <- sqrt(diag(vcovHC(mc40a)))

mc40b <- lm(propelig ~ splitcongswing40
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc40b <- sqrt(diag(vcovHC(mc40b)))

mc40c <- lm(propelig ~ pickupcongswing40
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 42),])
se.mc40c <- sqrt(diag(vcovHC(mc40c)))

mc42a <- lm(propelig ~ swingcong42
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 43),])
se.mc42a <- sqrt(diag(vcovHC(mc42a)))

mc42b <- lm(propelig ~ splitcongswing42
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 43),])
se.mc42b <- sqrt(diag(vcovHC(mc42b)))

mc42c <- lm(propelig ~ pickupcongswing42
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat[which(dat$year == 43),])
se.mc42c <- sqrt(diag(vcovHC(mc42c)))

# Output as tidyframe
tidycong <- rbind(tidy(mcavg) %>% 
                     filter(term == "swingcongavg") %>% 
                     mutate(model = "Swing = 50-55% (All Years)"),
                   tidy(mc36a) %>% 
                     filter(term == "swingcong36") %>% 
                     mutate(model = "Swing = 50-55% (1936)"),
                   tidy(mc36b) %>% 
                     filter(term == "splitcongswing36") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1936)"),
                   tidy(mc36c) %>% 
                     filter(term == "pickupcongswing36") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1936)"),
                   tidy(mc38a) %>% 
                     filter(term == "swingcong38") %>% 
                     mutate(model = "Swing = 50-55% (1938)"),
                   tidy(mc38b) %>% 
                     filter(term == "splitcongswing38") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1938)"),
                   tidy(mc38c) %>% 
                     filter(term == "pickupcongswing38") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1938)"),
                   tidy(mc40a) %>% 
                     filter(term == "swingcong40") %>% 
                     mutate(model = "Swing = 50-55% (1940)"),
                   tidy(mc40b) %>% 
                     filter(term == "splitcongswing40") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1940)"),
                   tidy(mc40c) %>% 
                     filter(term == "pickupcongswing40") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1940)"),
                   tidy(mc42a) %>% 
                     filter(term == "swingcong42") %>% 
                     mutate(model = "Swing = 50-55% (1942)"),
                   tidy(mc42b) %>% 
                     filter(term == "splitcongswing42") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1942)"),
                   tidy(mc42c) %>% 
                     filter(term == "pickupcongswing42") %>% 
                     mutate(model = "Swing = 47.5-52.5% (1942)"))
tidycong$std.error <- c(as.numeric(se.mcavg[2]),
                         as.numeric(se.mc36a[2]),
                         as.numeric(se.mc36b[2]),
                         as.numeric(se.mc36c[2]),
                         as.numeric(se.mc38a[2]),
                         as.numeric(se.mc38b[2]),
                         as.numeric(se.mc38c[2]),
                         as.numeric(se.mc40a[2]),
                         as.numeric(se.mc40b[2]),
                         as.numeric(se.mc40c[2]),
                         as.numeric(se.mc42a[2]),
                         as.numeric(se.mc42b[2]),
                         as.numeric(se.mc42c[2]))
tidycong$term <- c("Swing = 50-55% (All Years)", "Swing = 50-55% (1936)",
                    "Swing = 47.5-52.5% (1936)", "Swing = 45-50% (1936)",
                    "Swing = 50-55% (1938)", "Swing = 47.5-52.5% (1938)", 
                    "Swing = 45-50% (1938)", "Swing = 50-55% (1940)",
                    "Swing = 47.5-52.5% (1940)", "Swing = 45-50% (1940)",
                    "Swing = 50-55% (1942)", "Swing = 47.5-52.5% (1942)", 
                    "Swing = 45-50% (1942)")

# And plot
congress <- dwplot(tidycong, dodge_size = 0,
                   dot_args = list(size = 1.25),
                   whisker_args = list(size = 1),
                   vline = geom_vline(xintercept = 0, 
                                      colour = "grey60", 
                                      linetype = 2)) +
  theme_bw() + xlab("Coefficient Estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "none") + 
  scale_colour_grey(start = .25, end = .35)
# print #
png('congress.png', res = 120)
print(congress)
dev.off()


##### Table 2 -- Regional Results #####
mreg1 <- plm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + as.factor(year),
             data = dat[which(dat$year >= 41),],
             effect = "time", index = c("year"),
             model = c("random"))
se.mreg1 <- coeftest(mreg1, 
                     vcov=vcovHC(mreg1, type = "sss", cluster = "time"))[,2]

mreg2 <- plm(propelig ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + as.factor(region)
             + as.factor(year),
             data = dat[which(dat$year >= 41),],
             effect = "time", index = c("year"),
             model = c("random"))
se.mreg2 <- coeftest(mreg2, 
                     vcov=vcovHC(mreg2, type = "sss", cluster = "time"))[,2]

mreg42 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 42),])
se.mreg42 <- sqrt(diag(vcovHC(mreg42)))

mreg43 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 43),])
se.mreg43 <- sqrt(diag(vcovHC(mreg43)))

mreg44 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 44),])
se.mreg44 <- sqrt(diag(vcovHC(mreg44)))

mreg45 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 45),])
se.mreg45 <- sqrt(diag(vcovHC(mreg45)))


## Output to stargazer ##
stargazer(mreg1, mreg2, mreg42, mreg43, mreg44, mreg45,
          se = list(se.mreg1, se.mreg2, se.mreg42, se.mreg43, 
                    se.mreg44, se.mreg45),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Proportion Enlisted in County"),
          title = ("Impact of county presidential vote share on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County", "New England", 
                               "Confederacy", "County GOP Vote Share",
                               "Pct. Urban", "Total Population (1000s)",
                               "Crop Value", "Pct. Male Eligible", "Pct. Black",
                               "Farms/Population",
                               "Mid Atlantic", "North East Central", 
                               "West North Central", "South Atlantic", 
                               "East South Central", "West South Central",
                               "Mountain", "Pacific",
                               "1942", "1943", "1944", "1945",
                               "Intercept"))


##### Table 3 -- Predicting 1944 presidential election #####
swingm1 <- lm(demswing4044 ~ propelig + swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop
              + as.factor(year), 
              data = dat[which(dat$year >= 41 & dat$year <= 43),])
se.swingm1 <- sqrt(diag(vcovHC(swingm1)))

swingm2 <- lm(demswing4044 ~ propelig * swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop
              + as.factor(year), 
              data = dat[which(dat$year >= 41 & dat$year <= 43),])
se.swingm2 <- sqrt(diag(vcovHC(swingm2)))

swingm3 <- lm(demswing4044 ~ propelig + swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop, 
              data = dat[which(dat$year == 43),])
se.swingm3 <- sqrt(diag(vcovHC(swingm3)))

# and output
stargazer(swingm1, swingm2, swingm3,
          se = list(se.swingm1, se.swingm2, se.swingm3),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Democratic Vote Share Swing (1944-1940"),
          title = c("Associations between enlistment rates and vote share swing towards the Democratic presidential candidate, 1944 election. Robust standard errors in parentheses."),
          covariate.labels = c("Enlistment Rate", "Swing County",
                               "Pct. Urban", "Population (1000s)",
                               "Avg. Crop Value", "Pct Male (21+)",
                               "Pct. Black", "Pct. Farms",
                               "1942", "1943",
                               "Enlistment Rate * Swing County",
                               "Intercept"))